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ABSTRACT 

A new method for the inclusion of ionizing radiation from uniform radiation fields into 
3D Smoothed Particle Hydrodynamics (sphi) simulations is presented. We calculate 
the optical depth for the Lyman continuum radiation from the source towards the 
SPHI particles by ray-tracing integration. The time-dependent ionization rate equa- 
tion is then solved locally for the particles within the ionizing radiation field. Using 
test calculations, we explore the numerical behaviour of the code with respect to the 
implementation of the time-dependent ionization rate equation. We also test the cou- 
pling of the heating caused by the ionization to the hydrodynamical part of the sphi 
code. 

Key words: Methods: numerical - hydrodynamics - radiative transfer - Hn regions 



1 INTRODUCTION 



ZEUS (5tone fc Norman 1992) or nirvana (Ziegier, Yorke 



Kaisig 1996[), or codes includ i ng effects of IR and UV 



Smoot hed Particle Hydrodynamics (SPHj has become a nu- radiation (lYorke fc Kaisig 1995]; jonnhalter, Preibisch & 



merical method widely used for addressing problems related 
to fluid hows m astrophysics. Due to its Lagrangian nature 
it is especially well suited for applications involving varia- 
tions by many orders of magnitude in density. Examples for 
this type of applications are simulations of the collapse of 
molecular clouds and the formation of a stellar cluster, as 



Yorke 1993; Richling fc Yorke 1996 



performed by Klessen, Burkert & Bate (1998). A compar- 
ison between grid based methods and SPH was performed 
by Burker t, Bate & Bodenheimer (1996) and Bate & Burk- 



In contrast, the addition of physical processes to SPH 
codes is just at its beginnings. Extensions achieved so far 
are sophisticated equations of state (e.g. Vanhala et al. 1998) 
and self-gravity. Some efforts were made to make SPH faster 
and more accurate. The introduction of tree a lgorithms 
( |Barnes fc Hut 1989| |Press 198c] ; |Benz et al. 199C| ) and the 



ert (1997). They applied both methods to the numerically 
demanding problem ot gravitational collapse and tragmcn- 
tation of a slightly perturbed rotating cloud core with an 
r~ 2 density profile. Both m ethods yielded the same qual- 
itative results. Bate (1998) performed the first calculation 
which followed the collapse of a molecular cloud core in 3 
dimensions down to a protostellar object in hydrodynam- 
ical equilibrium, thus spanning 17 (!) orders of magnitude 
in density. Other applications include accretion pro c esses in 
massive circu mbinary disks ( Bonnell fc Bate 1994 ; Bate & 
Bonnell 1997), the collapse of clo ud cores induced by shock 
waves (Vanhala fc Cameron 1998 ) or colliding clumps (Bhat- 
tal et al. 199 j), he precession of accre tion disks in binary 



use of GRAvity PipE (grape), a hardware devi ce for fast 
computation of the gravita tional N-body forces ( Umemura 



et al. 1993; Steinmetz 1996), helped reducing the numerical 



effort for the gravitational force calculation and the deter- 
minat ion o f the nearest neighbours for each particle. Inut- 
suka (1995) presented a Godunov-like solver for the Eulerian 
equations in SPH thus enhancing the numerical treatment 
of sho cks. The introduction of gravitat i onal periodic bound- 
aries (Hernquist, Bouchet & Suto 1991; Klessen 1997) allows 



systems ( |Larwood fc Papaloizou 1997), the dynamical be- 
haviour of massive protostellar disks ([Nelson et al. 1996 ) or 



the treatment of fragmentation and turbulence in molecu- 
lar clouds without global collapse. The timestep problem 
which arises during isothermal collapse calculations at high 
densities is circumvented by the formation of sink particles, 
which substitute the innermost parts of the collapsing clump 
by one particle and accumulate the infalling mass and mo- 
menta (Bate, Bonnell & Price 1995). 



the formation of large scale st ructure and galaxies in the 
early universe (Steinmetz 1996). 



A variety of physical processes are at work in the in- 
terstellar medium, like magnetic fields, radiation or ther- 
mal conductivity, necessitating their inclusion into numeri- 
cal codes. This has already been achieved to a large extent in 
grid based methods like the magneto-hydrodynamics codes 



The strength of SPH lies in its Lagrangian nature, which 
makes it especially attractive for problems involving gravi- 
tational collapse a nd st ar formation. Applications like e.g. 
by Klessen et al. (1997), which deal with the collapse and 



fragmentation of molecular clouds, neglect the feedback pro- 
cesses of newly born stars which act on their parental cloud 
through stellar winds, outflows and ionization. This simpli- 
fication may be justified as long as the simulations deal with 
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collapse on timescales smaller than ~ 1 Myr, on which single 



/ Si l) a v dv 



and bi nary stars or T T au-like clusters are formed (Efrc- 
mov & Elmegreen 1998). The case is different for larger 



k = riH • a — riH 



timescales, on which OB subgroups and associations are 
formed. Neglecting feedback in these cases can lead to un- 
physical results, like a star formation efficiency of 100 per 
cent, since in the purely isothermal case all material will 
sooner or later be accreted onto the evolving protostel- 
lar cores. This is in strong contradiction to observations, 
which estimate a global star formation efficie ncy for or- 
dinary mo lecular clouds of order 10 per cent ( Wilking & 
Lada 1 )8E ) . Another possible effect of feedback is the in- 



duction of star formation due to the compression of cloud 
material by shock waves and ionization fronts. 

In this paper we discuss the implementation of the ef- 
fects of ionizing UV radiation by massive stars into SPH cal- 
culations as a first step in order to perform collapse cal- 
culations on scales where OB-stars are formed in a more 
realistic way. This will in future applications allow us to as- 
sess questions like: How does the process of ionization by 
massive stars change the stellar initial mass function? What 
are the implications for the star formation efficiency? Can 
star formation be induced by ionization, and if yes, what 
are the time scales and the parameter space, for which in- 
duced star formation can be expected? These questions will 
be discussed in subsequent papers. 



2 PHYSICAL PROBLEM 

We incorporate the effects of ionizing radiation from hot 
stellar photospheres into SPH by dividing the problem into 
three major substeps: 

(i) calculation of the UV radiation field by solving 
the time-independent, non-relativistic equation of radiative 
transfer, 

(ii) determination of the ionization and recombination 
rates from the local radiation field, density and ionization 
fraction, 

(iii) advancing the ionization state of the particles in time 
by solving the time-dependent ionization rate equation. 



9 (i) 

°tot 



(2) 



where er„ denotes the ionization cross section of hydrogen in 
the ground state and riH the particle density of the H atoms. 

The role of dust in Hn regio ns and its effect on ionizing 
radiation is still very uncertain ( Feldt et al. 1998 ). If dust 
is present, it will partially absorb UV photons, heat up and 
reemit the energy in the IR regime. Its first order effect can 
be included easily under the assumption of a homogeneous 
distribution of the dust in the Hn region. The corresponding 
contribution to the optical depth can be incorporated by 
adding the dust absorption coefficient at the Lyman border 
Kd to the absorption coefficient in Eq. |l|. depends on 
the dust model used and is regularity determined using Mie 
theory for grains with given distributions in size and shape. 
In this paper, we set «a to zero throughout. 

We also neglect the diffuse field of Lyman continuum 
photons, which are being produced by recombinations of 
electrons into the ground level and which themselves pos- 
sess sufficient energy for ionizing other H atoms. A thorough 
treatment of this radiation can only be achieved by detailed 
radiatio n tra nsfer calculations as proposed e.g. by Yorke & 
Kaisig (1995). Instead we use the assumption of the validity 
of the 'on the spot' approximation as follows: due to the fact 
that the spectrum of the Lyman recombination photons as 
well as the the ionization cross section is strongly peaked at 
the Lyman border, a small amount of H atoms in the ion- 
ized region is sufficient to make the medium optically thick 
for the Lyman recombination photons. This leads to the ab- 
sorption of these photons in the ultimate vicinity of their 
creation sites. As the creation of one photon is related to 
the creation of one H atom, its absorption leads to the de- 
struction of one H atom. Thus the net effect of these photons 
on the local ionization structure is zero. 

This assumption breaks down in regions next to OB 
stars, where due to the high UV flux the density of H atoms 
is not sufficient to make the medium optically thick to Ly- 
man continuum photons. Next to ionization fronts, where 
the density of H atoms is much higher, the 'on the spot' as- 
sumption is neverthele ss a g ood approximation. On further 
details refer to Yorke fll988). 



2.1 Calculation of the UV radiation field 

Given a planar infall of ionizing photons from a distant 
source onto the border of the volume of interest with a flux 
Jo Lyman continuum photons per time and square area, the 
resulting photon flux inside this volume is given by 

J(s) = Jo • exp (— t (s)) , 

where r(s) is the optical depth for the ionizing photons along 
the line of sight parallel to the infall direction of the pho- 
tons, and s is the distance from the border of the integration 
volume along the line of sight: 

t(s)= f [R(s')+K d (s)]ds'. (1) 
Jo 

We neglect the effect of 'photon hardening', i.e. the 
stronger absorption of weaker photons, and use an 'effec- 
tive' absorption coefficient R, the mean of k„ over frequency, 
weighted by the spectrum of the source S„: 



2.2 Ionization and recombination rates 

The ionization rate in the medium is given by the sinks of 
the UV radiation field, since every ionization leads to the 
absorption of one UV photon: 

X = n H aJ = -V ■ J, (3) 

where J = Je s is the flux vector in the direction e s of the 
line of sight. 

The recombination rate can be estimated as : 

1Z = n^as = n 2 x 2 aB, (4) 

with n being the particle density of H atoms and protons 
together, n e the particle density of free electrons, x = n c /n 
the ionization fraction and ob the effective recombination 
coefficient under assumption of validity of the 'on the spot' 
approximation. The recombination coefficient a is given as 
the sum over the individual recombination coefficients a n , 
where the electron ends up in the atomic level n: 
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(•») 



Under the assumption of the 'on the spot' approximation 
recombinations into the ground level do not lead to any net 
effect and thus a\ can be neglected in Eq. [|. The resulting 
net recombination rate which is used in Eq. ^ is commonly 
called tt B aft er the nomenclature introduced by Baker & 
Menzel fll96S|): 



qb = 



£< 



2.3 Ionization rate equation 

Knowing the ionization and recombination rates, 1 and 1Z, 
the ionization fraction can be calculated from the ionization 
rate equation. The time dependency of the ionization frac- 
tion in the frame comoving with the corresponding particle, 
i.e. its Lagrangian formulation, is given by: 

dn c 
"dT 



1 — TZ. 



(6) 



2.4 Modeling the source 



Since the spectral distribution of the UV radiation emitted 
by the photospheres of intermediate to high mass stars is 
very uncertain, we assume a black radiator with an effective 
temperature T*. 



3 NUMERICAL TREATMENT 

We developed two different methods for the numerical treat- 
ment of time dependent ionization in the SPH calculations. 
Both have in common the method of finding paths from 
the ionizing source to the particles, along which the optical 
depth for the Lyman continuum photons can be calculated. 
They differ in the way the ionization rate is determined given 
the radiation field. Method A uses the SPH formalism to 
calculate the divergence of the radiation field in Eq. ^. In 
method B we adopt a different approach also used in grid 
methods, where we derive the ionization rate from the dif- 
ference in the numbers of photons entering and leaving a 
particle. 



3.1 Finding the evaluation points on the path 
towards the source 

First, we specify the position, the rate of ionizing photons 
Stat and a (from Eq. 0) of the source. 

For each particle i we now proceed in the following way 
(see Fig. |l|): Given the list of nearest neighbours of particle 
i, which has to be determined anyway for the SPH formalism, 
we look for the particle j in the list, closest to the line of sight 
defined by the smallest angle O between the line connecting 
the particles i and j and the line of sight. We choose the angle 
between, not the distance from, the line of sight, since we are 
interested in controlling the error in the direction towards 
the source. This is not garanteed by the latter criterion. 

We store this particle in a list and determine the evalua- 
tion point Sj as the projected particle position on the line of 



sight. To determine the next evaluation point Sk even closer 
to the source we now repeat this method using the neighbor 
list of particle j and so forth until we reach the source. 

3.2 Calculating the optical depth and ionization 
rate for the particles 

3.2.1 Method A: SPH formalism method 

Now the path from the source to particle i is known, and 
the integration of Eq. (Q) can be discretized by using the 
evaluation points Si. The value for tih can be estimated by 
using the SPH smoothing formalism: 



tih (r) = nudW (r — n) 



(7) 



where the sum runs over the particle corresponding to the 
evaluation point and its nearest neighbours. W is the weight 
factor for each neighboring particle provided by the smooth- 
ing kernel. We calculate the optical depth along the line of 
sight by applying the Trapezian Formula, until we reach 
particle i: 



Tk+l 



1~k + -O" (Sk+l — Sk) (jlH,k+l + RH,k) ■ 



with Sk being the position of the evaluation point on the 
line of sight. Note that this treatment neglects the effects of 
scattering of the ionizing photons by recombination or dust. 

The distance between two successive evaluation points 
is smaller or equal to the local smoothing length, which de- 
termines the largest distance of the particles included in the 
nearest neighbour list as well as the spatial resolution. This 
guarantees that the line of sight integration of Eq. ([!]) is 
discretized into a reasonable amount of substeps, consistent 
with the resolution given by the underlying particle distri- 
bution. 

The flux of ionizing photons at the position of particle i 
into the direction of photon propagation e s is then given by: 

Ji = Jo • e s ■ exp (— t (s)) . 

With the ionizing flux known at the particle positions, 
the nabla operator in Eq. (|^) can be calculated by the SPH 
formalism. It is given for each particle i as the sum over its 
neighbours: 



^ Pi 



Now we are able to solve Eq. ^, which we write as: 



— — = I; — TliX: qb- 

at 



(8) 



(9) 



The time scale for the establishment of ionization equilib- 
rium is given by 1/ (nets ) , which is regularly much shorter 
than the dynamical and gravitational timescales we are in- 
terested in. In order to avoid small timesteps arising from 
the usage of explicit methods, we use an implicit scheme. 
The first order discretization of Eq. ^| over a time interval 
At is given by: 



n+1 n a , 

X? — x\ + At ■ 



rn+l 



n+1 n+1 \ 

n s x i a B ) , 



(10) 



where the indices n and n+1 denote the values at the begin- 
ning and the end of the actual timestep At, respectively. We 
already know all the values on the right hand side from ad- 
vancing the particles by the SPH formalism, except the value 
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Figure 1. Illustration of the path finding procedure. Each plus sign represents a particle. The circle segments symbolize the radius of the 
volume filled with particles in the nearest neighbour list of the corresponding particle. The particle with the smallest angle between 
the line of sight from SO to S5 and the line connecting them with the target are used for the determination of the evaluation points for 
the integration along the line of sight (small circles) . 



for Therefore a fully consistent implicit treatment is 

not feasible. We use the following guess for this value: 

■exp(-< +1 aaf +1 (l-xf +1 )) 



Xf 



;2f 



1 - exp (-nf +1 aaf +1 (1 - x! 1 )) 



(11) 



In this equation, we assign an effective radius ai to each 
particle i proportional to the mean particle separation, given 
by di = (Mj/pi) 1 ' 3 . This is the estimate of the size of a region 
with the particle mass Mi and density pi. The factor with the 
exponentials on the right hand side accounts for the effect of 
higher absorption and hence ionization rate with decreasing 
ionization fraction. 

We must use the effective radius ai in Eq. [If] instead of 
the smoothing length h, since the method works analogous 
to implementations in grid codes. In contrast to the SPH 
formalism, each particle now represents a volume of total 
mass Mi and density pi, in which ionizing radiation enters 
on one side and leaves on the opposite side. The size of this 
volume is given by ai as defined above. It is proportional to 
the particle spacing. 

In contrast, hi differs from the mean particle separation 
as it is defined by the condition that there is a fixed number 
of neighbors -/V ne igh of mass M in the sphere with radius 2 hi 
and is thus given as 



hi = 



3N ncigh M 
32irp 



1/3 



It depends on iVneigh and can therefore not be used instead 
of ctj in Eq. |l| 

One consequence of the discretization of the ionization 
rate equation is that the solution in ionized regions tends 
to oscillate around the equilibrium value. In order to avoid 
small timesteps arising from this, we set the ionization frac- 
tion x of particles with an x > 0.95 to the equilibrium value 
ie, which is defined by dn c /dt = in Eq. ^: 

da; 1 dn B 2 

37 = 7T = f(l - xe)J - nx E a B = 0. 

at n at 

With k — aJ/(naB) follows that 



IE 



Af + 4k 



1/2 



- k 



This method works well in absolutely smooth, noise free 
particle distributions. However, if one wishes to initially dis- 
tribute the particles randomly in space, one runs into prob- 
lems. The sum in Eq. ^ is very sensitive to noisy particle 
distributions. Eventually the noise can be so high, that the 
error of the sum introduced by noise reaches the order of the 
sum itself. The ionization rate then locally drops below zero 
for some particles, which can only be avoided by smoothing 
the ionization rate spatially over several smoothing lengths. 
The result is poor resolution. We circumvent this problem 
in method B. 



3.2.2 Method B: grid based method 

In this case, a different method is used to discretize the cal- 
culation of the optical depth. We determine the positions of 
the evaluation points i along the line of sight as described 
in Sect, ^t] and calculate the hydrogen density iiH.i at these 
positions using Eq. |. The path is then divided into pieces 
with length Asi = (si+i — sj_i)/2, assuming a constant hy- 
drogen density nu.i along each interval. The optical depth 
for one piece can then be approximated by 

Ati = (TTlH.iAsi. 

These contributions to the optical depth are summed up 
until we reach the position located one effective radius Oi 
before the position of particle k. A first order approximation 
for the ionization rate is now given by 



Jo 



2a k n k 



exp (Yk-a) (1 - exp (-Art)) , 



where Tk_ a = ^2 Ari denotes the optical depth one effective 
radius before the particle's position and Ar k = 2aknH,k^" 
the optical depth across the particle. 

With the ionization rate derived above we solve the ion- 
ization rate equation as described for case A. One can easily 
show that Eqs. ( p"o| ) and (|Tl]) now give the exact implicit first 
order discretization for Eq. (9). The solution now approaches 
the equilibrium value xe in the ionized regions without the 
instabilities mentioned in method A. It is not necessary to 
set x artificially to xe- 
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Figure 2. Mean number of evaluation steps / per path depend- 
ing on tolerance angle © t ol an d number of particles TV. Upper 
panel: / depending on to i for different JV. Note how / drops 
with increasing to i- Lower panel: log-log-plot of / depending on 
TV for different © to i with the scaling law / <x TV 1 / 3 overplotted 
as a solid line. For to i > 0°, I becomes independent of TV for 
large TV. 



Method A seems to be the more consistent method since 
it uses the SPH formalism for the calculation of I. This is 
the reason why it is also discussed in this paper. Neverthe- 
less we prefer method B due to its robustness against noisy 
particle distributions and higher consistency concerning the 
integration scheme and have applied it to a couple of test 
cases. 



the source to Sj . Thus no integration is needed for this part 
of the path. One only has to make sure that Tj is already 
known, i.e. that the line of sight integration for particle j has 
been performed earlier. In this case, the average number of 
evaluation points I per line of sight only depends on O to i for 
large TV. As shown in Fig. ^, I becomes constant for large 
TV and decreases with increasing Otoi- As soon as / becomes 
independent of TV the total computational effort for all lines 
of sight together scales tx TV. 

We demonstrate the effects of using the tolerance angle 
on the accuracy of the ionization rate calculation in Fig. |3[ 
Histograms are plotted for the errors in T and r for calcula- 
tions with Otoi = 0.5° , 1° , 2° and 90° compared to O to i = 0° . 
As the particle distribution we chose the evolved state of a 
numerical simulation which studies the compression and col- 
lapse of a dense clump within the UV field of an OB associ- 
ation using 200 000 particles. The results of this calc ulation 
will be presented elsewhere (Kessel & Burkert 1999). Note 
that Otoi = 90° represents the worst case, since the toler- 
ance angle criterion now is fulfilled for every particle with 
minimal O per search through the nearest neighbour list. 

The particles which are most affected by the tolerance 
angle criterion lie next to the borders of shadows cast by op- 
tically thick regions, since here the path for the integration 
along the line of sight may be bent through the optically 
thick region, thus decreasing the ionizing flux artificially. In 
the opposite case, the path may lead around the opaque 
region, increasing the ionizing flux at the position of a par- 
ticle in the shadow. These extreme cases lead to the tail in 
the error histograms in Fig. |j. Applying the tolerance angle 
criterion thus numerically blurs shadows. 

The mean errors in r are 1.3 per cent for O to i = 0.5°, 
2.2 per cent for O to i = 1.0°, 3.4 per cent for O to i = 2.0° and 
11.2 per cent for O to i = 90°. The correspnding mean errors 
in X are 2.8, 4.1, 5.7 and 13.3 per cent, respectively. For 
the remaining test cases presented in this paper the choice 
of Otoi has no effect, since they deal with one-dimensional 
problems, in which the optical depth is only a function of 
distance from the source. Applying the tolerance angle cri- 
terion only shifts the evaluation points from the lines of 
sight in directions perpendicular to these, along which there 
is no change in the optical depth. Indeed even the choice 
Otoi = 90° gives the same results in the one-dimensional 
test cases as for O to i = 0°. Thus the errors introduced by 
the angle criterion must be checked with problems in which 
this symmetry is broken and shadows are present, as the one 
mentioned above. 



3.3 Computational effort 

If the procedure explained above is used, the computational 
effort for the line of sight integration scales approximately 
as TV 4//3 , since the integration has to be done for each of the 
TV particles, and the average number of evaluation points on 
each line of sight scales as TV 1 / 3 . 

We can reduce the exponent from 4/3 to 1 by intro- 
ducing a 'tolerance angle' Otoi- Suppose we determine the 
particles along the line of sight as expalined . As soon as O 
for a particle j along the line of sight towards the source is 
smaller than O to i we stop our search here. The optical depth 
of this particle Tj is then used as an estimate of the opti- 
cal depth along the remaining part of the line of sight from 



3.4 Smoothing the ionization front 

For reasons of noise reduction we smooth the ionization 
front, which is not resolvable by the SPH representation, over 
a distance of the order of one local smoothing length. Na- 
ture provides a simple way for doing this. The width of the 
ionization zone is of the order of one photon mean free path 
length, 

d= (cT-nn)-\ (12) 

where a is the net absorption cross section for ionizing pho- 
tons as defined in Eq. [| 

Since we cannot resolve the ionization region anyway, 
we are free to adjust a in a way that the width of the ion- 
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Figure 3. Histograms of the relative errors in r and X for different 
© to i in a three-dimensional test case. 

ization region given by Eq. [L2| is equal to a constant factor 
C < 1 times the local smoothing length h, but never larger 
than the value a given by Eq. ^: 

a — min \a, (nn • C • /i)" 1 ] . 

Test calculations have shown that a good value is C = 0.1. 
It has proven to sufficiently reduce numerical noise intro- 
duced into the ionization structure by noise in the particle 
distribution and at the same time to keep the resolution of 
ionization fronts better than the resolution of the SPH for- 
malism in order not to worsen the overall resolution. Note 
that, when "smoothing" the ionization front over 0.1 times 
the smoothing length, the noise reducing effect is not caused 
by the spatial smoothing, since it is ten times smaller than 
the SPH smoothing. It rather results from a larger number 
of time steps needed to ionize a particle in the front from an 
ionization fraction of x = to a; ~ 1. This gives the neigh- 
bouring particles the opportunity to react to the changed 
state in a smoother way. 

3.5 Heating effect 

We assume that heating and cooling effects lead to an equi- 
librium temperature of 10 000 K in the ionized gas pene- 
trated by ionizing radiation. The cross sections for elastic 
electron-electron and electron-proton scattering are of the 
order 10 _13 cm 2 . Together with a mean velocity of the elec- 



trons of the order of 600 km s -1 the thermalization timescale 
for the energies of the ejected electrons is far less than a year 
for densities of 1 particle cm -3 , which is many orders of 
magnitude smaller than the dynamical timescale. Thermal- 
ization thus occurs quasi instantaneously. This process runs 
even more rapidly for higher densities. Thus we are allowed 
to treat the gas behind the ionization front as thermalized. 
We set the internal energy to: 

e — x ■ eioooo + (1 — x) • e co id, 

with eioooo being the internal energy corresponding to a tem- 
perature of 10 000 K for ionized hydrogen, and e co id to the 
internal energy for the 10 K cold, neutral gas. Note that this 
method does not properly treat recombination zones, since 
in this case one needs the correct inclusion of the heating and 
cooling processes in order to achieve the correct gas temper- 
atures, sound velocities and pressures. Also, the equilibrium 
temperature in Hll-regions can vary by 20 per cent from this 
value. These deviations can also only be taken into account 
by proper treatment of heating and cooling. 



4 TESTS OF THE NUMERICAL TREATMENT 

Although being of one-dimensional nature, the following test 
problems were performed fully in three dimensions. 

4.1 Test 1: Ionization of a slab with constant 
density 

With this problem we test the implementation of the time- 
dependent ionization rate equation by ionizing a slab of Hi 
gas of constant density n with ionizing radiation falling per- 
pendicular onto one of the boundary surfaces. With hydro- 
dynamics switched off, we let the ionization front traverse 
the slab with a constant velocity V{. To achieve this, we have 
to vary the infalling photon flux with time. It is given by 

J(t) = Jf -+- Jt = nvf + n 2 aBVft, 

where the first term on the right hand side is the flux which 
provides the photons being absorbed in the ionization front. 
The second, time-dependent term equals the loss of photons 
on their way through the slab until they reach the front. 

For the initial setup we place a number N of parti- 
cles randomly into a slab with length-to-height and length- 
to-width ratios of 10. Subsequently we let the particle dis- 
tribution relax by evolving it isothermally within the slab, 
adding a damping term to the force law. This is necessary to 
diminish the numerical noise which was introduced by the 
random distribution. We now have an ensemble of the par- 
ticles which does not possess any privileged directions and 
which represents a gas of constant density and temperature. 
We use this distribution as our starting configuration. From 
now on we keep the particles fixed in space and switch off 
hydrodynamics. 

The test was performed for a total number of N = 
2 000,16 000 and 128 000 particles. Since the spatial reso- 
lution for SPH calculations scales as A r_1/ ' 3 (with number of 
neighbours N na i s h per particle fixed), this yields an increase 
of linear resolution of a factor of two from one simulation 
to the simulation with next higher resolution. The results of 
these tests are shown in Fig. W. 
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Figure 4. Relative error in ionized mass vs. time between calcu- 
lations and theoretical result for test 1. Plus signs: 2000 particles, 
stars: 16000 particles, diamonds: 128000 particles. Time in units 
of time needed for the ionization front to cross the slab of length 
^slab with a propagation velocity V{. 



The mean relative errors between the theoretical result 
and the calculations decrease linearly with increasing reso- 
lution, consistent with our first order discretization of both 
the line of sight integration and the time dependent ion- 
ization equation. The error also decreases with time as the 
representation of the ionization front gets thinner and thin- 
ner compared to the already ionized region. The spread in 
the errors for N = 16000 and N = 128000 results from 
the fact that in these cases the numerical solution oscillates 
around the theoretical solution, sometimes being larger than 
the latter, sometimes smaller. 



4.2 Test 2: Ionization of a slab with density 
gradient 

We proceed as in test 1, with the difference that we choose 
a slab with a constant density gradient in the direction of 
photon propagation. We choose the time dependence of J 
such that the ionization front should travel through the gas 
with constant V{. J is given by: 



T ... / 2 dn \ 

J(t) = n v{ + I Q B rio + jT u f I 



dn 2 2 

V{t + a>BTlo—Vft + 

ax 

fdn\ 2 3 3 



where no denotes the density at the surface where the radi- 
ation penetrates the slab and 4 s is the density gradient. 

In Fig. |E] we plot the ionized mass for the theoretical 
solution and the numerical simulations against time. The 
numerical results converge against the theoretical solution 
with increasing resolution. The deviations at t > 0.9 are 
caused by the ionization front reaching the rear boundary 
of the slab. 

Note that the version of SPHI used in this paper is not 
able to follow ionization fronts exactly which travel faster 
than one local smoothing length per time step. This must 
be taken into account during the timestep determination. 
In applications with fast ionization fronts (typically R-type 
fronts in the early phases of the evolution of Hn-regions) this 
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Figure 5. Ionized mass vs. time for test 2. Solid line: theoret- 
ical solution. Plus signs: N=2 000. Stars: N=16 000. Diamonds: 
N=128 000. 
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Figure 6. Density profile for test case 3 and different resolutions 
of the SPHI calculation. Unit of the x-axis normalized to the po- 
sition of the shock front. Ionizing radiation infall from the left. A 
shock wave traveling to the right into the undisturbed medium 
with no = 10 cm -3 and T co i^ = 100K sweeps up a dense shell of 
post-shock material, which is separated from the thin, hot, ion- 
ized material by an ionization front. Solid line: analytical result. 
Ratio of the thickness of the swept-up layer to the current local 
smoothing length for different resolutions: dotted 6, dashed 12, 
dash-dot 20. Corresponding times in code units: 0.34, 0.70, 1.0. 



criterion can lead to very small timesteps and thus to a high 
amount of CPU time needed. A version which circumvents 
this problem is being developed. 



4.3 Test 3: Coupling of ionization and 
hydrodynamics 

For this t est w e adopt the problem mentioned by Lefloch & 
Lazareff (1994). A box filled with atomic hydrogen of par- 
s and temperature T co id = 100 K 
with the photon flux in- 



ticle density no = 10 cm~ 
is exposed to ionizing radiation, 
creasing from zero linearly with time with a rate d$/dt = 
5.07 ■ 10 -2 cm _2 s -2 . There exists an analytical solution to 
this problem, which is self similar in the sense that physical 
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Figure 7. Velocity profile for test case 3 and different resolutions. 
Unit of the x-axis normalized to the position of the shock front. 
Ionizing radiation infall from the left. Solid line: analytical result. 
Ratio of the thickness of the swept-up layer to the current local 
smoothing length for different resolutions: dotted 6, dashed 12, 
dash-dot 20. Corresponding times in code units: 0.34, 0.70, 0.93, 
1.0. 



values at position x measured in the direction of the pho- 
ton flow at time t are only functions of x/t. This means: 
the structure is stretched with time. The convergence of 
the code towards the correct solution with increasing resolu- 
tion can be tested in one calculation, since for all appearing 
structures the ratio between structure sizes and smoothing 
lengths increases linearly with time. 

The resulting structure is the following: an isothermal 
shock is driven into the neutral medium, sweeping up a dense 
layer of material. This is followed by an ionization front 
which leaves the ionized material in quasi-static equilib- 
rium (see Figs. ^[Q). Using the parameter A = a~ 1 (d$/dt), 
Lefloch & Lazareff ( L994 ) find the following analytical solu- 
tion: 



A = nfVi 
> i 

n A 2 



A 



no 



Vi = 



Act 



\n J 



where rii, no and n c denote the particle densities of the ion- 
ized gas, the undisturbed neutral gas and the gas in the 
compressed layer, respectively, and Vi and V s the velocities 
of the ionization front and the shock front, respectively. 

We adopt a B = 2.7 ■ HT 13 cm 3 s" 1 from Lefloch & 



Lazareff (1994) in order to directly compare the results of 
the SPHI code to those of their grid-based method using a 
piecewise l inear scheme for the advection terms proposed by 
Van Leer (1979). The resolution of 192 grid cells along the 



slab of their calculations, from which they derived their re- 
sults, is comparable with the one used in our high resolution 
case. We use the same method as described in Sect. 4.S to 



produce the initial conditions. No gas is allowed to enter or 
leave the surface. 

Table ^ lists the result of this comparison. The SPHI 
calculation slightly underestimates V s and Vi, as is also ob- 



Table 1 . Comparison of analytical and numerical results for test 
case 3. 



analytical SPHI 



Lefloch e.a. (1994) 



rii (cm 3 ) 
n c (cm -3 ) 
V B (km s" 1 ) 
V (km s" 1 ) 



0.756 
1.59 ■ 10 2 
3.71 
3.48 



0.75 ± 0.05 
(1.55 ± 0.05) 
3.67 ±0.05 
3.43 ± 0.05 



10- 



0.748 
1.69 ■ 10 2 
3.51 
3.36 
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Figure 8. Density of the layer for test case 3. Unit of the x-axis 
normalized to the position of the shock front. Ionizing radiation 
infall from the left. Solid line: analytical result. Ratio of the thick- 
ness of the swept-up layer to the current local smoothing length 
for different resolutions: dotted 6, dashed 12, dash-dot 18, dash- 
dot-dot-dot 20. Corresponding times in code units: 0.34, 0.70, 
0.93, 1.0. 



served for the grid code. The errors of order 5 per c ent ar e 
comparable to those achieved by Lefloch & Lazareff (1994). 



In the early phases, i.e. low resolution, the poor treat- 
ment of the ionization front leads to irregularities in the ion- 
ized region and thus produces sound waves travelling back 
and forth between the boundary to the left and the ioniza- 
tion front (Figs. ^, |^), which decrease in power as time in- 
creases, i.e. at higher resolution. With increasing resolution, 
i.e. increasing ratio of layer thickness to smoothing length, 
the representation of the dense layer and the shock front 
improves (Figs. B, 0). 



5 SUMMARY 

The method presented in this paper allows the treatment 
of the dynamical effects of ionizing radiation in SPH calcu- 
lations. Thus the study of astrophysical problems arising 
from ionization, like the impact of ionizing radiation from 
newly born stars onto the evolution of their parental molec- 
ular clouds or the more consistent treatment of heating by 
OB associations in galaxy dynamics calculations are now 
feasible for the first time with SPHI in 3 dimensions. We 
demonstrate that the code is able to treat time- dependent 
ionization, the related heating effects and hydrodynamics 
correctly. Our first applications, detailed calculations of pho- 
toionization induced collapse in molecular clouds and results 
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Figure 9. Velocity profile for test case 3. Unit of the x-axis nor- 
malized to the position of the shock front. Ionizing radiation infall 
from the left. Solid line: analytical result. Ratio of the thickness 
of the swept-up layer to the current local smoothing length for 
different resolutions: dotted 6, dashed 12, dash-dot 18, dash-dot- 
dot-dot 20. Corresponding times in code units: 0.34, 0.70, 0.93, 
1.0. 

obtained from them, will be presented in a subsequent pa- 
per. 

To allow the correct treatment of recombination zones, 
one has to include the effects of time dependent heating and 
cooling processes by ionization and recombination, emission 
of forbidden lines and thermal radiation from dust. Another 
important aspect which was neglected here is the effect of the 
diffuse Lyman continuum recombination field. It can lead to 
the penetration of regions shielded from the direct ionizing 
radiation by the ionization front, which is e.g. seen in cal- 



culations o f photoevaporating prot ostellar disks (Yorke & 
Welz lfgl |Richling fc Yorke 199c]). An implementation of 



Efremov Y. N., Elmegreen B. G., 1998, MNRAS, 299, 588 
Feldt M., Stecklum B., Henning Th., Hayward T. L., Lehmann 

Th., Klein R., 1998, A&A, 339, 759 
Hernquist L., Bouchet F. R., Suto, Y., 1991, ApJS, 75, 231 
Inutsuka S.-L, 1995, IAUS, 174, 373 
Kessel, O., Burkert, A., 1999, in preparation 
Klessen R. S., 1997, MNRAS, 292, 11 

Klessen R. S., Burkert A., Bate M. R., 1998, ApJL, 501, 205 
Larwood J. D., Papaloizou J. C. B., 1997, MNRAS, 285, 288 
Lefloch B., Lazareff B., 1994, A&A, 289, 559 
Nelson A. F., Benz W., Adams F. C, Arnett D., 1998, ApJ, 502, 
342 

Press W. H., 1986, in Springer Lecture Notes in Physics Vol. 267, 
The Use of Supercomputers in Stellar Dynamics, Springer, 
Berlin 

Richling S., Yorke H. W., 1998, A&A, 340, 508 

Sonnhalter C, Preibisch T., Yorke H. W., 1995, A&A, 299, 545 

Steinmetz M., 1996, MNRAS, 278, 1005 

Stone J. M., Norman M. L., 1992, ApJS, 80, 791 

Umemura M., Fukushige T., Makino J., Ebisuzaki T., Sugimoto 

D., Turner E. L., Loeb A., 1993, PASJ, 45, 311 
Vanhala, H. A. T., Cameron, A. G. W., 1998, ApJ, 508, 291 
Van Leer, B., 1979, J. Comp. Phys., 32, 101 

Wilking, B. A., Lada, C. J., 1985, in Black D. C, Matthews M. 
S., eds., Protostars & Planets II, The University of Arizona 
Press, p. 297 

Yorke, H. W., 1988, in Chmielewski Y., Lanz T., eds., Radiation in 
Moving Gaseous Media, Advanced Course of the Swiss Society 
of Astrophysics and Astronomy vol. 18, pp. 242—246. Geneva 
Observatory, CH-1290 Sauverny-Versoix, Switzerland. 

Yorke H. W., Kaisig M., 1995, Comp. Phys. Comm., 89, 29 

Yorke H. W., Welz A., 1996, A&A, 315, 555 

Ziegler U., Yorke H. W., Kaisig M., 1996, A&A, 305, 114 



these processes into our SPHI code is planned in the future. 
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